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Abstract 

A  numerical  method  for  a  finite  difference  approach  has  been  established  for  the  analysis 
and  control  of  the  fluid  field  behavior  of  flow  past  a  cylinder.  The  discretization  of  the  2D 
Navier-Stokes  equations  is  done  over  a  staggered  grid,  the  convective  terms  in  the  momentum 
equations  are  handled  using  a  mixture  of  central  differences  and  donor-cell  discretization,  and 
the  Poisson  equation  for  the  pressure  is  solved  through  the  successive  overrelaxation  (SOR) 
method. 

We  also  study  some  open-loop  and  closed-loop  control  of  the  flow  field  by  rotating  the 
cylinder.  For  the  open-loop  control  design,  we  mainly  make  use  of  the  energy  method,  and 
for  the  feedback  design  both  the  energy  method  and  the  phase  lock-in  method  are  applied. 

Keywords:  Finite  Difference;  Flow  Control;  Vortex  Shedding. 

Extended  Abstract 


1  Introduction 

Numerical  computation  is  an  indispensable  tool  for  control.  Computation  is  needed  both  for  off¬ 
line  functions  such  as  simulation,  analysis,  and  synthesis  of  control  systems  as  well  as  for  on-line 

*  Paper  is  submitted  to  the  33rd  AIAA  Fluid  Dynamics  Conference  and  Exhibit,  Orlando,  Florida  23  -  26  Jun 
2003 

Research  supported  in  part  by  the  Summer  Faculty  Fellowship  Program  from  Air  Force  Office  of  Scientific 
Research,  2002. 


functions  associated  with  control  system  implementations  in  embedded  processors.  Simulation 
is  used  to  extend  our  capacity  to  reproduce  and  forecast  physical  processes  on  the  computer. 
Expensive  experiments  are  increasingly  being  replaced  by  computer  simulations.  Moreover, 
simulation  enables  the  examination  of  processes  that  cannot  be  experimentally  tested. 

Strouhal  was  the  first  to  study  the  periodic  features  produced  by  the  movement  of  a  cylinder 
body  in  air  at  low  Reynolds  number  [18].  The  understanding  of  the  spatiotemporal  dynamics  of 
oscillatory  flows  has  proceeded  by  successively  considering  the  linearized  version  of  Navier-Stokes 
equations  by,  among  others,  Joseph  [17],  Monkewitz,  Huerre  and  Chomaz  [21],  and  Yudovich 
[29].  Since  the  Navier-Stokes  equations  linearized  about  a  steady-state  flow  is  non-self-adjoint, 
the  investigation  is  a  difficult  problem,  and  for  the  solution  of  the  majority  of  questions  it  is 
necessary  to  use  numerical  methods  [29] . 

In  this  paper,  we  use  the  finite  difference  approach  to  investigate  wake  responses  under  open- 
loop  and  closed-loop  control.  Simulations  for  controlling  the  wake  flow  by  rotating  the  cylinder 
are  provided. 


2  Problem  Description 


The  flow  behavior  in  general  can  be  described  more  or  less  as  follows.  Let  u(t)  denote  the 
velocity  field  of  the  wake  flow  (here  u(f)  is  the  symbolic  notation  for  a  function  of  the  space 
variable  £  =  (£i , £2. £3))-  When  Re  is  small  (Re  <  40),  the  velocity  field  is  stable  and  attracts 
all  of  the  orbits.  In  this  case  the  flow  is  fully  laminar.  When  40  <  Re  <  50,  steady  vortices  form 
behind  the  cylinder.  Typically,  for  t  ->  00,  u(f)  will  converge  to  a  stationary  solution 

u(t)  — >  us,  as  t  — >  00, 

the  limit  solution  depending  on  the  initial  condition  u0.  In  this  case  the  flow  is  still  steady. 

When  Re  is  grater  than  approximately  50,  a  Hopf  bifurcation  occurs,  i.e. 

u(t)  -  ip(t)  — >  0  as  t  — >  00,  (2.1) 

where  ip  is  a  time-periodic  solution  of  the  system.  After  the  Hopf  bifurcation  has  occurred, 
the  flow  never  becomes  stationary  again.  In  fact  it  is  time  periodic.  In  this  case  von  Karman 
vortices  move  in  the  flow  direction  and  new  vorticies  appear  on  the  downstream  face  of  cylinder 
in  a  seemingly  time-periodic  manner,  namely,  vortex  shedding. 

The  system  of  equations  governing  the  conservation  of  momentum  for  fluid  flow  is  known 
as  Navier-Stokes  equations.  When  combined  with  conservation  of  mass,  they  describes  the 
general  evolution  of  fluid  flow.  For  a  viscous  incompressible  fluid,  the  equations  are 

( du  \ 

P  +  (U  •  V)uJ  -  pAu  +  Vp  =  f  (2.2) 

div  u  =  0,  (2.3) 

where  p  >  0  is  the  constant  density  of  the  fluid,  v  >  0  is  the  kinematic  viscosity,  and  /  represents 
volume  forces  that  are  applied  to  the  fluid. 


2 


Alternatively,  equation  (2.2)  can  be  considered  as  the  nondimensionalized  form  of  the  Navier- 
Stokes  equations  in  which  case  u,  p,  f  are  nondimensionalized  quantities  and  the  kinematic 
viscosity  p/p  :=  v  is  replaced  by  Re*1,  Re  representing  the  Reynolds  number 


where  U  and  L  are  the  typical  velocity  and  length  used  for  the  nondimensionalization. 

In  this  paper,  we  only  consider  the  two-dimensional  case  and  rewrite  equations  (2.2)  in 
component  form.  Let  u  =  (u,v),  p  =  pp,  and  f  =  p(fx,  fy).  Then  the  equations  in  component 
form  read  as  follows: 

momentum  equations: 

du  d(u 2)  d(uv)  1  ( d2u  d2«\  _ 

dt  Ox  ~  dy  Re  ^  dx"2,  dy2  J  x 

dv  d(uv)  d(v2)  1  / d2v  d2v\  _ 

dt  dx  dy  Re  ^ckc2  dy2  J  y’ 


(2.4) 

(2.5) 


continuity  equation: 


du  dv 
dx  dy 


(2.6) 


Let  Cl  be  the  region  exterior  to  a  circular  disk  with  radius  r.  Then  the  open-loop  or  closed-loop 
control  is  acting  on  the  circumference: 


(u,v)  =  g(t)  on  dCl. 


(2.7) 


3  Finite  Difference  Approach 


Our  numerical  calculation  is  made  over  a  collection  of  uniformly  space  discrete  grid  points.  When 
central  differences  are  used  for  the  incompressible  Navier-Stokes  equations,  pressure  oscillations 
can  possibly  occur  if  all  three  unknown  functions  u,  v  and  p,  are  evaluated  at  the  same  grid 
points  [2].  There  are  two  ways  suggested  for  fixing  the  problem.  One  is  to  utilize  the  upwind 
differences,  and  the  other  is  to  use  a  staggered  grid.  The  advantage  of  using  the  staggered  grid 
is  that  a  central  difference  approach  can  still  be  maintained.  A  staggered  grid  is  illustrated  in 
Fig.l.  The  pressures  are  calculated  at  the  cell  center,  the  horizontal  velocities  u  are  evaluated  in 
the  midpoints  of  the  vertical  cell  edges,  and  the  vertical  velocities  are  computed  in  the  midpoints 
of  horizontal  cell  edges.  Cell  ( i,j )  occupies  the  spatial  region  [(i  -  l)Sx,i8x]  x  \(j  -  l)5y,jSy], 
and  the  corresponding  index  (i,  j)  is  assigned  to  the  pressure  at  the  cell  center  as  well  as  to  the 
u- velocity  at  the  right  edge  and  the  v- velocity  at  the  upper  edge  of  this  cell.  The  key  feature  is 
that  the  pressures  and  velocities  are  evaluated  at  different  grid  points  so  that  the  non-physical 
oscillations  due  to  the  numerical  computation  can  be  avoided. 
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Staggered  Grid 


Figure  1:  Stagged  grid  cells. 

3.1  Discretization  of  the  spatial  derivatives 
Let  the  rectangular  region  be 

ft  =  [0,o]  x  [0,6] 

on  which  we  introduce  a  grid.  The  grid  is  divided  into  imax  cells  of  equal  size  in  the  ^-direction 
and  jmax  cells  in  the  ^-direction.  We  denote 

Sx  =  — -Q-  ,  8y  =  -  -  . 

lmax  Jmax 

Then  using  the  central  difference  scheme,  we  have 


d2u  | 

ui+l  ,j  “ 

2  Uij  ~f  — 1 J 

(3.8) 

dx 2  1  ij 

(<Sx)2 

d2u  | 

+ 1  “ 

2  MiJ  + 

(3.9) 

dy 2  It j' 

to 

d2v  | 

vi±l,j  ~~ 

2^t,j  lj 

(3.10) 

dx 2 

(<Sx)2 

d2v  | 

vi,j+ 1  — 

2wtj  +  Vij-i 

(3.11) 

dy2  1  i,j 

(<h/)2 

dp\ 

Pi+l,j  “ 

PiJ 

(3.12) 

dx  1  ij 

Sx 

dp\ 

dyUj 

PiJ+ 1  ~ 
8y 

PiJ 

(3.13) 

4 


3.2  Discretization  of  convective  terms 


Considering  that  the  convective  terms  in  the  momentum  equations  can  become  dominant  at  high 
Reynolds  numbers,  following  [11],  we  use  a  mixture  of  the  central  differences  and  the  donor-cell 
discretization.  In  each  cell  we  set 


d(u2) 
dx  i,j 


J_  +  Ui+lj'j2  _  +  Uij 

,  1  (  I  ui,j  Uj+ljl  (Uitj  —  Ui+ lj)  +  Uij  I  (Ui-Ij  —  Uij ) 

+7fe  l - 2 - 2  2  2 


d(uv) | 

1 

/  ivi,j  “5"  vi+l,j)  ui,j+ 1) 

(Ui,j- 1  +  Vj+lJ-l)  (Uij—i  +  Uij)\ 

dy  1  i,j 

8y 

-L'-v. 

V  2  2 

1  f\vij  +  Ui+ij  (Uij  —  Uij+ 1) 

2  2  J 

\Uij-l  +  Uij. l,j-l|  {ui,j- 1  —  ui,j) 

5y  V  2  2  2  2 


9(uu) | 

(ui-ij  4-  Uj-ij+i)  (v<-ij  +  vi,j)\ 

dx  1  ij 

V  2  2 

2  2  ) 

,  1  ^  |Wjj  +  Ui  j4-i  |  (Ujj  Ui+ij)  \Ui— ij  +  Uj—ij+il  (Vi—ij  Uij) 

+7fe  l - 2  2  2  2 


^(f2) |  =  i  ({vu  +  vij+i\2  ( 2'ii-i  +  r\ 

dy  I* j  <5y  \\  2  /  V  2  /  / 

,  1  (  \Ui,i  +  ^.i+l I  Kj  -  «i J+l)  Kj-1  +  vi,j I  Kj-1  -  vij)\ 

+1Jy  \ - 2  2  2  2  )  ’ 

where  7  is  a  parameter  between  0  and  1.  During  the  simulation,  the  7  is  chosen  to  be  0.9. 


3.3  Iteration  for  solving  the  Poission  equation  for  pressure 


We  use  a  forward  difference  for  the  time  discretization.  Thus  the  momentum  equations  become: 


Let 


u(n+l)  =  u(n)  +  5t 


1 

Re 


f  d2u  d2u\ 

y&c2  dy2  J 


y(n+l)  =  y(n)  +  St 


1  /  d2v 

Re  yc>£2  dy2  J 


d(u2) 

dx 

d(uv ) 
dx 


d(uv) 

dy 


d(f) 

dy 


+  fy 


F(n)  =  u(n)  +  St 


1 

Re 


( d2u  d2u\ 
ycfa2  dy2) 


G(n)  =  v{n)  +  St 


1 

Re 


< 


d2v  d2v\ 

dx2  dy2  J 


d(u2)  d(uv)  i  t 

~Yx  dy~  +  U 

d(uv)  d(v2)  t 

~d~x  W  +  iy 


(3.14) 

(3.15) 


(3.16) 

(3.17) 
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We  then  set 


u(n+1>  =  F(n)  -  St 


Qp(n+ !) 


t/(n+1)  =  G{n)  -  St 


dx 
dp<'n+1') 
dy 


(3.18) 

(3.19) 


Making  use  of  the  continuity  equation,  we  have  the  Possion  equation  for  the  pressure  p(n+U  at 
time  tn+i : 


02p(n+l)  ^  02p(n+l)  j  C>G(n)N\ 

dx 2  ^  dy 2  St  y  dx  ^  dy  J 

Using  central  difference  for  the  Laplacian  operator,  we  obtain 

Jn+l)  _  2v{n+1)  +  o(n+1)  o(n+1)  -  2n(n+1)  +  Jn+1) 

p»+i,j  ,  Pij+i  zPt,j  +Pi,j- 1 

(dx)2  +  (dy)2 

p(») 


1  /  PW  _  /'tW  /'M") 

_  1  (  fjj  ^i-1.7  f  Gi,j  ~  Gt,j-1 


dt 


dx 


Sy 


(3.20) 

(3.21) 

(3.22) 


An  iterative  method,  called  successive  overrelaxation  (SOR),  is  applied  to  solve  the  above  dis¬ 
cretized  equation: 


Pfj1  =  Pij  +  u 


- - p'f. 

+  _2_ 


,  ( W  ^  (to 


(rjt  .  4.  r)^+l  nit  I  /  r>(n)  Z^C71)  r*(n)  /~i(n)  > 

Pt+ij+Pt-i,j  ,  Pij+i  +  Pij-1  1  (  i’ij  -  Pi-1, j  ,  -Gjj.j 
(dx)2  (dy)2  dt\^  dx  +  Ty 

where  u  €  [0, 2]  is  a  parameter.  The  residual  is  defined  to  be 


-«  _  & U  -  +  rf-,  j  ,  P&+,  -  2pfj  +  K&_.  1  /iff  -  ,  G*?  -  (*>_,  \ 

iJ  (^)2  +  (Sy)2  Tt[  Si  +  Ty - j 

The  iteration  is  terminated  when  either  a  specified  number  of  iterations  has  been  reached  or  the 
norm  of  the  residual  is  smaller  than  a  specified  tolerance.  In  the  simulation,  u  =  1.7. 


3.4  Stability  issue 


One  of  common  methods  to  ensure  stability  of  the  numerical  algorithm  and  prevent  the  numerical 
solution  from  generating  non-physical  oscillations  is  to  set  the  time  step  to  meet  three  conditions: 

5x  .  Sy 

<  1 - 7,  St  <  ;  —  ~ .  (3.23) 


ju  Re  (  1  1  V1 

<  2  {(Sx)2  + (Sy)2)  ’  5t 


\'Umax\ 

We  follow  reference  [25]  and  use  an  adaptive  stepsize  control: 

Sx 


\Vn 


St  :=  rmin 


(Re/  1  1  V1  6x 

[2  \(Sx)2  +  (Sy)2)  *|W| 


Sy 


I  Umax  I  \vmax  |  J 

where  the  factor  r  G  [0, 1]  is  a  safety  factor.  In  the  simulation,  we  choose  r  =  0.5. 


(3.24) 
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4  Simulation  with  Open-loop  Control 

The  simulations  are  based  on  the  solution  of  the  following  Navier-Stokes  equations: 


^  =  -Vp-( u  V)u  +  4~V2u  in  Q,  (4-25) 

ot  rie 

div  u  =  0  in  fi,  (4-26) 

u  =  s(x2,—xi)/r  on<9f l,  (4-27) 

u  =  Uoo,  onTi  (4.28) 

u  =  0  onr2,r3,  (4.29) 


where  s  is  the  angular  speed  (possibly  time-dependent)  of  the  cylinder,  and  =  (u^,  0)  is 
the  free  stream  velocity. 


5  Simulation  with  Closed-loop  Control 


It  is  well-known  that  the  von  Karman  vortex  street  induces  a  strong  unsteady  transverse  periodic 
force  on  the  cylinder  which  can  lead  to  cylinder  oscillations  with  large  amplitudes.  Such  oscil¬ 
lations  can  cause  structural  damage  and  structure  failure  in  applications.  A  means  of  altering 
the  wake  structure  through  feedback  control  in  order  to  reduce  or  eliminate  such  vibrations  is  of 
practical  importance.  One  approach  is  to  introduce  feedback  into  the  system  with  one  or  more 
sensors  and  actuators  so  that  the  combined  system  is  stable.  Here  we  study  a  simple  approach. 

We  assume  that  we  can  measure  the  upstream  velocity  away  from  the  cylinder.  In  our 
simulations,  we  use  the  average  of  the  speeds  of  three  cells  (which  are  about  1.5  units  ahead  of  the 
cylinder)  to  serve  as  known  information.  Then  we  place  five  sensors  evenly  on  the  downstream 
semicircle  of  the  cylinder.  Let  us  denote  these  five  sensing  velocities  as  ( ).  Let  the  active 
control  be 


g(f)  =  Gsin(27r/t)(y,  —  x)/r  on  dfl  (5.30) 

where  G  is  the  feedback  gain,  r  is  the  radius  of  the  cylinder,  and  /  is  a  selected  frequency. 
Here  the  frequency  /  can  be  chosen  such  that  it  can  reduce  the  vibration  of  the  cylinder.  For 
example,  /  can  be  picked  to  be  larger  than  the  shedding  frequency  to  avoid  the  occurrence  of 
resonance.  The  feedback  gain  G  is  set  to  be  the  following.  Let  T  =  1/f  and 

A (tn)  :=  E  VM*n)  -  Ui(tn  -  !T))2  +  M*n)  -  «i(*n  ~  ^))2.  (5.31) 

i 

Then  at  t  =  0  we  set  G  =  U^.  When  the  feedback  is  turned  on,  after  one  period  T  we  set 

G(t)  =  G(tn)  +  A(tn)  for  tn  <  t  <  tn+i  :=  T  +  tn  (5.32) 

if  the  free  stream  velocity  remains  the  same.  Otherwise  reset  G  to  be  the  free  stream  velocity. 
The  purpose  is  to  let  the  feedback  increase  the  amplitude  automatically  until  the  flow  is  locked 
in  a  periodic  solution  with  the  period  T.  In  this  case  the  cylinder  will  be  oscillated  at  a  desirable 
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frequency.  The  numerical  simulations  of  the  closed- loop  design  are  carried  out.  The  simulations 
are  done  in  multiple  phases.  [The  selected  frequency  (or  forced  frequency)  is  /  =  0.25.]  The 
natural  frequency  (without  the  control)  is  about  0.21.  In  phase  1,  [/«>  =  1  and  Re  =  100,  and  in 
phase  2,  U^o  =  5  and  Re  =  500.  In  phase  3,  Ur^  =  1  and  Re  —  100,  which  are  the  same  values 
as  in  phase  1. 


Phase  Lock-in  Feedback  Control 


Figure  2:  Phase  Lock-in  Feedback 


6  Outline  of  Simulation  Results 

A  sequence  of  numerical  simulations  have  been  obtained.  The  simulations  consists  of  the  flow 
behaviors  (1)  without  control;  (2)  with  open-loop  control;  (3)  with  closed-loop  control.  We  here 
only  provide  some  results  of  the  simulations  (Fig.3-6). 
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Streaklines  at  time  50.0486 


Figure  3:  Streakline  in  a  vortex  shedding  period  without  control 


Vortictty  vibes  *t  lime  30.1 


Figure  4:  Vorticity  distribution  at  Re  =  100  and  time  «  30  when  the  cylinder  is  rotated  in  a 
constant  speed  =1 
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